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We have developed a time propagation scheme for the Kadanoff-Baym equations for general 
inhomogeneous systems. These equations describe the time evolution of the nonequilibrium Green 
function for interacting many-body systems in the presence of time-dependent external fields. The 
external fields are treated nonperturbatively whereas the many-body interactions are incorporated 
perturbatively using ^-derivable self-energy approximations that guarantee the satisfaction of the 
macroscopic conservation laws of the system. These approximations are discussed in detail for the 
time- dependent Hartree-Fock, the second Born and the GW approximation. 



I 

O 

o 



> 
o 

o 
o 



X 



PACS numbers: 31.15.xm, 31.15.-p 



I. INTRODUCTION 



The recent developments in the field of molecular elec- 
tronics have emphasized the need for further development 
of theoretical methods that allow for a systematic study 
of dynamical processes like relaxation and decoherence 
at the nanoscale. Understanding these processes is of ut- 
most importance for making progress in molecular elec- 
tronics, whose ultimate goal is to minimize the size and 
maximize the speed of integrated devices To study 
these phenomena, theoretical methods must allow for the 
possibility to study the ultrafast transient dynamics 0, Q 
up to the picosecond [1, Q and femtosecond timescale, 
while including Coulomb interactions, without violating 
basic conservation laws such as the continuity equation 
0. A theoretical framework that incorporates these 
features is the nonequilibrium Green function approach 
based on the real-time propagation of the Kadanoff- 
Baym (KB) equations 0, 11, S tlO, [H [H, [11 [3 ■ This 
method allows for systematic inclusion of electron in- 
teractions while providing results in agreement with the 
macroscopic conservation laws of the system [1, @ . In two 
recent Letters 0, [HI we applied the KB equations to in- 
vestigate the short time dynamics of atoms and molecules 
in time-dependent external fields, as well as the transport 
dynamics of double quantum dot devices. It is the aim 
of this paper to describe in detail the underlying method 
that was only briefly described in those Letters. This in- 
cludes both a description of the theory as well as the time- 
propagation algorithm. We further generalize the equi- 
librium method, described in two recent papers [Tsl. [l6j. 
to the nonequilibrium domain. We also extend earlier 
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work on the time-propagation method of the KB equa- 
tions for homogeneous systems [13, [3 to the case of 
inhomogeneous systems. In the inhomogeneous case we 
can not take advantage of Fourier transform techniques 
anymore. The KB equations become time-dependent ma- 
trix equations instead, in which the matrices are indexed 
by basis function indices. The time-stepping algorithm 
has to take into account the special double-time structure 
of the equations which are furthermore nonlinear, inho- 
mogeneous and non-Hermitian. Therefore, several stan- 
dard time-propagation methods can not be used. Our ap- 
proach is different from the one presented in Refs. (l7l.[l8j 
by incorporating correlated initial states and the memory 
thereof, which is described in terms of Green functions 
with mixed real and imaginary time arguments. To sim- 
plify the time-stepping procedure, we make use of several 
symmetry relations of the Green function. 
The paper is divided as follows: in section II we present 
the KB equations and their symmetry properties. In sec- 
tion III we discuss the conserving self-energy approxi- 
mations that we use, and in section IV we present the 
time-propagation method that we developed for systems 
described within a general basis set representation. Fi- 
nally in section V we present a summary and conclusions. 



II. THEORY 

We consider a many-body system that is initially in 
equilibrium at a temperature T and with a chemical po- 
tential /i. At an initial time the system is exposed to a 
time-dependent external field. This external field can, for 
instance, be a bias voltage in a quantum transport case, 
or a laser pulse. The field forces the system out of equi- 
librium and we aim to describe the time-evolution of this 
nonequilibrium state. In second quantization the time- 
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dependent Hamiltonian of the system reads (throughout 
this paper we use atomic units h = m = e = 1) 



H{t)^ j dx V>t(x)/i(x,t)?A(x) + 

+ \ j j rfxidx2?/'^(xi)'0t(x2)i;(ri,r2)?/'(x2)-0(xiX,l) 

where x = (r, a) denotes the space- and spin coordinates. 
The two-body interaction will, in general, be a Coulombic 
repulsion of the form w(ri,r2) = l/|ri — Y-i\. The one- 
body part of the Hamiltonian is 



/i(x,<) = -iv2 + u;(x,i) 



(2) 



where w^Xji) is a time-dependent external potential. 
The chemical potential fj, of the initial equilibrium sys- 
tem is absorbed in the one-body part of the Hamiltonian. 
The expectation value of an operator O, for a system ini- 
tially in thermodynamic equilibrium {t < to), is given 
by 



(O) = Tr {pO}, 



(3) 



where p = e /'^o/Tr e is the statistical operator, 
Hq is the time-independent Hamiltonian that describes 
the system before the time-dependent field is applied and 
(3 — l/ksT is the inverse temperature. The trace here 
represents a summation over a complete set of states in 
Fock space [l9]. After the time-dependent external field 
is switched on at time to, the expectation value is given 
by 



(0(t)) = 



Tr{f7(to-z/3,to)Og(t)} 
Tr {[/(to -1/3, to)} ' 



(4) 



where Onit) = U{tQ,t)OU{t,tQ) is the opera- 
tor O in the Heisenberg picture and C/(t2,ti) = 
T[exp(— i J^^ dtH(t))], for t2 > ti, is the time-ordered 
evolution operator of the system. We further wrote 
exp(— /3iJo) = U(to — iP,to) as an evolution operator in 
imaginary time. If we read the time arguments in Eq.(l4|) 
from right to left we see that they follow a time-contour 
as displayed in Fig[Tl This contour is also known as the 
Keldysh contour [20l. [2l| . A more detailed inspection of 
Eq. (HI then shows that the expectation value can also be 
written as a contour-ordered product [S^, [H, [U, [2^ . 
The one-particle Green function is then defined as a 
countour-ordered product of a creation and an annihi- 
lation operator 



G(l,2) = -i(Tc[V'//(l)ViJ^(2)]), 



(5) 



where Tc denotes the time-ordering operator on the con- 
tour and where we used the compact notation 1 = (xi, ti) 
and 2 = (x2,t2). If we consider the Green function at 
time ti = to — iP and use the cyclic property of the trace, 
we find that G(xito - 2) = -G(xito,2) [23|. Hence, 
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FIG. 1: Keldysh contour. The depicted contour allows for 
the calculation of observables for times to < t < T. The 
initial Green function is calculated on the imaginary track 
[to, to — As we propagate the KB equations in time, for 
real times t > to, the turning point of the time-contour at 
t = T moves to the right along the real time axis. 



the Green function defined in Eq. ([5]) obeys the boundary 
conditions 

G(xito,2) = -G(xlto-^/3,2), (6) 
G(l,X2to) = -G(l,X2to-^/3). (7) 

The Green function satisfies the equation of motion 

[idt, - h{l)]G{l,2) ^ 5(1,2) + / d3E(l,3)G(3,2), (8) 

Jc 

as well as a corresponding adjoint equation In 
Eq.® the time-integration is carried out along the con- 
tour C. The self-energy E incorporates the effects of ex- 
change and correlation in many-particle systems and is a 
functional of the Green function that can be defined dia- 
grammatically 0, [l^ . The Green function can be written 
as 

G(l,2) = 0(t,t')G>(l,2) + 0(t',t)G<(l,2), (9) 

where is a step function generalized to arguments on the 
contour i.e. with 9(t,t') = 1 if t is later on the contour 
than t' and zero otherwise @. The greater and lesser 
components G^ and G^ respectively, have the explicit 
form 



G>(l,2) = -z(V'i/(l)V;l,(2)), 
G<(l,2) = z(7^1,(2)^ff(l)). 



(10) 
(11) 



When one of the arguments is on the vertical track of the 
contour, we adopt the notation [l^l 

Gl(l,X2,-^r2) = G<(l,X2,to-ir2), (12) 
Gr(xi,-iri,2) = G>(xi,to -iri,2). (13) 

Finally, for the case when both time arguments are on 
the imaginary track of the contour, we have the so-called 
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Matsubara Green function iG^ [T^ 

iG*^(xiri,X2T2) = G{xito - iTi,X2to - ira), (14) 

which is a well-known object from the equilibrium theory. 
The factor i in the definition of Eq. p^ is a convention 
which ensures that G*^ is a real function. The self-energy 
S has a similar general structure as the Green function 

E(l, 2) = 2) + 9{t, t')^> (1, 2) + 9{t', t)S< (1, 2). 

(15) 

The main difference with Eq. Q is the appearance of the 
term T,^^ which is proportional to a contour delta func- 
tion 5{ti,t2) in the time coordinates [9]. This term has 
the explicit form 



[G](l,2) = (5(ti,i2)S^^(xi,X2,ti), (16) 



where 



E^^(Xi , X2, i) = lG< (Xit, X2t)«(xi , X2) 

-i^(xi-X2) J dx3w(xi,X3)G<(x3t,X3i). (17) 

The structure of this self-energy is that of the Hartree- 
Fock (HF) approximation. However, in general we will 
evaluate this expression for Green functions G obtained 
beyond HF level (see section Using the form of 

the self-energy of Eq. (fl5] ) the contour integrations can 
be readily carried out [9, 20] and we find separate equa- 
tions for the different Green functions G^, G^ ^ and G^^. 
To display their temporal structure more clearly we sur- 
press the spatial indices of the Green functions and self- 
energies. Alternatively, these quantities may be regarded 
as matrices [7|. On the imaginary track of the contour 
we obtain 



[^dr - h]G''' (r) = 5{t) + / dfE''' (r - f)G*' (f), (18) 

Jo 

where the Green function and the self-energy are func- 
tions of the time-differences only, i.e. iG^^ti — T2) = 
G{—iTi,—iT2) and zS*^(ti — T2) = S(— iti, — «T2), since 
the Hamiltonian is time-independent (and equal to Hq) 
on the imaginary track. Equation p^ . which determines 
the Green function of the equilibrium system, has been 
treated in detail in references [H, [lB|- For the other 
Green functions we obtain 

idtG^{t,t') = h"''{t)G^{t,t')+lf{t,t'), (19) 



-idt'G^{t,t') ^ G^{t,t')h"'' {t') + I^{t,t'), (20) 
«9tGl {t, -ir) = h"^{t)G^ {t, -It) + /I (<, -ir), (21) 
~idtG^{~iT, t) = G^{-iT, t)h"^{t) + I^-ir, t), (22) 

where h"^{t) = h{t) + I]^^(i) and E^^(i) is given by 
Eq.([T7]). The retarded and advanced functions for G and 
I] are defined according to 

i^«M(t,t') = ±ei±tTt')[F>{t,t') - F<{t,t')l (23) 



with F replaced by G and E respectively. The so-called 
collision terms 1$ and T have the form 



lfit,t' 



= f dt^'^{t,t)G^{t,t')+ f dtY.^{t,t)G^{t,t') 
Jo Jo 

1 r'^ 

+ - dfY){t,-if)G^{-if,t'), (24) 
* Jo 

= f dtG^{t,t)T,^{t,t')+ f dtG^it,t)^'^{t,t') 
Jo Jo 

1 

+ - / dfG\t,-if)j:H-if,t'), 

« Jo 



l](t,-iT)^ / dtJ:^{t,t)G\t,-iT) 



(25) 



+ / dfT,^t,~if)G^{f -t), (26) 
Jo 

I^{-iT,t)^ [ dtG^{-iT,t)^^{t,t) 

Jo 

+ / dfG^^ir ~f)j:^{~if,t). (27) 

These equations are readily derived using the conversion 
table of Ref. [1^ . From the symmetry relations 

G^{t,t')^ = -G^{t',t), (28) 

E^(t,t')^ = -S^(t',t), (29) 

it follows that we only need to calculate G'* (i, t') and 
E>(i,f) for t > t' and G<{t,t') and Y.<{t,t') for t < t' . 

These equations imply that I^2{^t^') — 0^- We 

further have 



G^{-lT,t)=G\t,-l{P-T))\ 

j:^{^iT,t) = Ei(t, -i(/3-T))t. 



(30) 
(31) 



The symmetry relations ((28)) and (|30p for the Green func- 
tion follow directly from its definition, whereas the sym- 
metry relations and (pij) for the self-energy follow 
from Eqs.(3.19) and (3.20) of Ref. fg"]. Another conse- 
quence of equations ([30)1 and ([3T|) is that I^{—iT,t) — 

[I^t,~i{f3-T))y , which means that in practice it is 
sufficient to calculate only and l\ Egs. lfT^ to 

are known as the Kadanoff-Baym equations @, Q ■ 
Once the Matsubara Green function G*^(t) is obtained 
from Eq. pS]) . the Green functions G^{x =^,] [) can be 
calculated by time propagation. Their initial conditions 
are 



G>(0,0) =^G*^(0+), 
G<(0,0) =iG*^(0-), 
Gl(0,-iT) =^G*^(-T), 
Gr(-zr,0) =^G*^(T). 



(32) 
(33) 
(34) 
(35) 



The KB equations, together with the initial conditions, 
completely determine the Green functions for all times 
once a choice for the self-energy has been made. The form 
of the self-energy will be the topic of the next section. 
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III. SELF-ENERGY APPROXIMATIONS 

In the applications of the KB equations it is possible to 
guarantee that the macroscopic conservation laws, such 
as those of particle, momentum and energy conservation, 
are obeyed. Baym Q has shown that this is the case 
whenever the self-energy is obtained from a functional 
$[G], such that 



E(l,2) = 



(5$ 



(36) 



Such approximations to the self-energy are called con- 
serving or ^-derivable approximations. Well-known con- 
serving^ approximations are the Hartree-Fock, the second 
Born [8], the GW [H], and the T-matrix [8] approxi- 
mation. In our work we implemented the first three of 
these. 

The second Born approximation - This approximation 
for the self-energy consists of the two diagrams to second 
order in the two-particle interaction [27| 

S(l,2) = S^^(1,2) + S(2)(l,2), (37) 

where is the HF part of the self energy of Eq. ^TB]) 

and E(2) = I](2a) + Y.i'^b) is the sum of the two terms 

E(2°) (1,2) = -i^G{l,2) j d3d4i;(l,3) 

xG(3,4)G(4,3)i;(4,2), (38) 



E(2'') (1,2)^^2 / d3d4G(l,3)i;(l,4)G(3,4) 



xG(4,2)t;(3,2), 



(39) 



where f(l,2) — ^(xi, X2)(5(ii, <2)- These terms are usu- 
ally referred to as the second-order direct and exchange 
terms. This approximation to the self-energy has been 
discussed in detail for the equilibrium case in Ref. p^ . 
For the nonequilibrium case we need to calculate the vari- 
ous components Yi^{x ] [). These are explicitly given 

by 

Y,(2a)S (1,2) = -i2G^(l,2) J d3rf4w(l,3) 

xG^^(3,4)G«^(4,3)i;(4,2), (40) 

j.(2a),ir (1^2) = -z^y d3d4Gir(l,2)i;(l,3) 

xGir(3,4)Gn(4,3)t;(4,2), (41) 
for the direct diagram, and 

Et^b).:? (1,2) = i'^ J d3d4G^(l,3)i;(l,4)G^(3,4) 

xG;§(4,2)t;(3,2), (42) 
E(26).l r (1, 2) = i^ J di dA GT r(l, 3)i;(l, 4)Gn (3, 4) 

xGir(4,2)i;(3,2), (43) 



for the second-order exchange diagram. These expres- 
sions follow immediately from Eos . (1551) and ([55]) with 
help of the conversion table of Ref. [20] . 
The GW approximation - In the GW approximation the 
exchange-correlation part of the self-energy is given as 
a product of the Green function G with a dynamically 
screened interaction W . The screened interaction 
W satisfies the equation 

IF(1,2) = i;(l,2) -I- J d3d4i;(l,3)P(3,4)M/^(4,2)(44) 

Here, v is the bare Coulomb interaction, and 

P(l,2)--iG(l,2)G(2,l), (45) 

is the irreducible polarization (26j . However, since the 
first term in Ea. (|44|) is singular in time (proportional to 
a delta function) it is convenient, for numerical purposes, 
to define its time-nonlocal part W = W — v [1^1 . From 
Eq.dMI) it follows that 



W{1,2)^ y d3d4w(l,3)P(3,4)i;(4,2) 

+ / d3rf4u(l,3)P(3, 4)11^(4, 2). (46) 



In terms of W, the self-energy has the form [20| 

E(l, 2) = E^^(l, 2) + iG(l, 2)H^(1, 2). (47) 

The part Ec = iGW represents the correlation part of 
the self-energy and has the components 

Ef(l, 2) = iG^(l, 2)1^^(2,1), (48) 
E]r(l,2) = iGTr(l,2)H^n(2,i). (49) 

From the fact that T^(l, 2) has the same symmetries 
as the contour-ordered density response function [2^ 
x(l,2) = — z(Tc[n/f (l)nH(2)]), where fi is the density 
operator, it follows that 

W^{2,1) = W^il,2) = -[H^^(2,l)]*, (50) 
(1, 2) = 1^^(2,1). (51) 

In the following, we will again surpress the spatial co- 
ordinates in order to display the temporal structure of 
the equations more clearly. From the symmetry relations 
pi)) . ([5T|) . ((i5|) and (gH), and the fact that we only need 
E>(t, t') for t > t' and E<(t, t') for t < t', it follows that 
we only need to calculate W'^ {t, —ir), and W^{t,t') for 
t < t' . The latter obey the equations: 

W<it, t') = vP<{t, t')v + vX<{t, t'), (52) 
(t, -ir) = vPT (<, -iT)v + vX^ {t, -ir), (53) 



where 



P<{t,t') = ~iG<{t,t')G>{t',t), (54) 
pl {t, -ir) = ~iG^ (t, -iT)GH~iT, t), (55) 
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and where the terms X"^ and X'^ are given by 

X<{t,t') = /* dtP<{t,l)W^{t,t') 
Jo 

+ [ diP^{t,i)W<it,t') 
Jo 

+ / dfP'^{t,~if)W^{-if,t'), (56) 

"'0 

X^t,~iT) = [ diP'^{t,i)W'^{t,~iT) 
Jo 

+ / dfP\t,-if)W^'\f -t), (57) 
Jo 

with the retarded and advanced quantities defined as in 
Eq. (|23l) . The initial conditions for and are 

W<iO,0) ^iW'^'iO-), (58) 
(0, -ir) = iW^\-T), (59) 

where iW^*^(T — T') — W(to—iT, tQ — iT') is the Matsubara 
interaction discussed in detail in Ref. n&i . 



IV. TIME-PROPAGATION OF THE 
KADANOFF-BAYM EQUATIONS 

In the following, we will describe the time-propagation 
method which we employed to solve the KB equations. 
This method can be applied to general Hamiltonians con- 
taining one- and two-body interactions, and is further in- 
dependent of the explicit form of the self-energy. 

The time-propagation method is applied to the KB 
equations in matrix form. This matrix form is obtained 
by expressing the Green function in terms of a set of ba- 
sis functions (/>i(x), which we choose to be Hartree-Fock 
orbitals 0, [ll, [S] 

G(xt,x't') = ^Gy(t,i')'/'^(x)0*(x'). (60) 

ij 

When Eq. ([60| is inserted in the expressions for the self- 
energy we obtain a basis set representation of the self- 
energy involving the matrices Gij{t,t') and the two- 
electron integrals which are given as integrals of orbital 
products with the two-body interaction v. All the quan- 
tities therefore become time-dependent matrices and all 
products are to be interpreted as matrix products. We 
will, however, surpress all matrix indices to display the 
temporal structure of the equations more clearly. Ex- 
plicit expressions of the matrix form of the second Born 
and GW self-energy are given in Refs.pTl. [T5l.[T6j. 
We start by discussing the time-propagation of G> and 
. Due to the symmetry relations Eq. ([28]) and ((29|) we 
only need to calculate G^{t,t') for t > t' and G^{t,t') 
for t <t'. From Eqs.dHD and ^ it then follows that 
G^ must be time-stepped in the first time-argument 




T T+A t 



FIG. 2: Time-stepping in the {t,t') -plane. G^{t,t') is calcu- 
lated for t > t' and G'"{t,t') is calculated for t < t' . 

and in the second one. We thus need to calcu- 
late G>{T + A,t') and G<{t,T + A) for a smaU time 
step A, from the knowledge of G^ (t, t') for t, t' < T. 
The symmetry relations (|28|) then immediately provide 
us with G>{t',T + A) and G<{T + A,t) as well. The 
time-stepping procedure is illustrated in Fig[2]that dis- 
plays the {t, t')-plane, in which at a given time T all the 
quantities inside the square with sides equal to T, are 
known. The time-step G<{t,T) ^ G<{t,T + A) corre- 
sponds to a shift of the upper side of the time square with 
A i.e. a shift from the solid to the dotted line in FiglS] 
Similarly the time-step G> (T, t') G> (T + A, t') corre- 
sponds to a shift of the righthand side of the time square 
with A. We further need to make a step G^{T,T) 
G<{T + A,T + A) along the time diagonal t = t'. The 
propagation of G^ {—ir, t) and G^ (t, —ir) requires a time- 
step in the real time coordinate t for fixed imaginary time 
points T. 

Note that the righthand sides of Egs. p^ to ([^ depend 
on the Green functions at the times T + A, which are not 
known at time T. We therefore carry out the time-step 
T ^ T + A twice. After taking the time step for the first 
time, we recalculate the righthand sides of Egs. fTO]) to 
(1^^ and repeat the time-step T ^ T + A using an aver- 
age of the old and new collision and HF terms. Since the 
term h^^{t) in Eas. (fT9| to (p2| can attain large values, it 
is favorable to eliminate this term from the time-stepping 
equations. For each time-step T ^ T + A we therefore 
absorb the term in a time-evolution operator of the form 

U{t) = e-*''""m*, (61) 

where h"^{T) = h{T + A/2) + S^^(T), where h is the 
one-body part of the Hamiltonian of Eq.®. The one- 
body Hamiltonian h{t) is explicitly known as a function 
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of time and can be evaluated at half the time-step. The 
term is only known at time T and will be recalcu- 
lated in the repeated time-step. In terms of the opera- 
tor U{t) of we define new Green function matrices 
g'^ix =^,10, as 

G^{hM) = U{h)g^{hM)U\t2), (62) 
Gl (ti, ~jT2) - U{h)g^ ih, ^iT2), (63) 
G^{-in,h) = i~in,h)uHh). (64) 



We can now transform Eqs. fTO)) to ([221) into equations 
for g^. For instance, g-^ satisfies the equation 

idtg>{t,t') = U^{t)ih"^it) - h"^)G>{t,t')Uit') 

+U\t)I>{t,t')U{t'). (65) 

Since h"^ « h"^{t) for times T < t < T 4- A, we can 
neglect for these times the first term on the right hand 
side of Eq. (|65)) . We then find 

- U{T + A)g> {T + A,t2)U\t2) = 

pT+A 

g>{T,h)+ / dtdtg>{t,h) 




e^''"'(*-^)|/i>(r,t2) 

^U{A)G>{TM) - V{A)I>{TM),. 
where V{IS) is defined as 



(66) 



(67) 



F(A) = {h^')-\l^e 



Similarly for G< , which is propagated using Eq. ([20]) , we 
find the equation 



G<(ti,T + A) 



G<(ti,T)t/t(A) 
-/2<(ti,T)0(A). 



For time-stepping along the time-diagonal we use 

idtG<{t,t)^[h"^{t),G<{t,t)] 

+I<(t,t)-I<{t,t), 



(68) 



(69) 



which follows directly from a combination of the equa- 
tions for G^ of Eqs. (fl9|) and ([20]) . The corresponding 
equation for g^{t,t) on the time diagonal then becomes 

tdtg< {t,t)^ f/t (t) [h"^ (t) - h"^ ,G<{t,t)] U{t) 

+U\t){I<{t,t)-I<{t,t))U{t). (70) 

From this equation we then obtain 

G<(T + A,r + A) = 
= [/(r + A).g<(T-F A,r-F A)f/1'(r + A) 
= U(A)G<{T,T)U\A) - 

.A 

iU{A) / dtU\t)h2U{t) C/t(A), (71) 



where we defined /12 = I^{T,T) - I<{T,T). By using 
the operator expansion 



e^Be-^ = B+[A,B] + [A, B]] 
+ \[a\[A, [A, 5]]] + ..., 



it follows that 



i / dtu\t)h2U{t) = X^c*^ 



where 



iA 

n + 1 



(72) 



(73) 



(74) 



and G(") = -iA/12. If we insert Eq.dTS]) into Eq.(l7Tl) we 
finally obtain 



G<(T-^A,T+A) = U{A) 



G<(r, T)-f ^g(") 



t/t(A) 
(75) 

We found that keeping terms for n < 3 only, yields suf- 
ficient accuracy. We now consider the time propagation 
for the mixed real and imaginary time Green functions. 
For g^ we have the equation 

idtgHt,-iT) ^ U{t)\h"^{t)~h"'")G'\ {t,-iT) 

+U{t)n\t,-iT). (76) 

This yields, similarly as in Eg. lpS)) and (p5)) 

Gl (T + A, -iT2) = C/(A)g1 (T, -iT2) 

-V{A)i\t,-zt2). {11) 

Finally, for G^ we have 

G^ {-iTi.T + A) = G^ {-iTi,T)U{A)^ 

-I^{-iTi,T)V{Ay. (78) 

The Eqs. JMl), JMl), (El), dZZl) and (UHl) form the basis 
of the time-stepping algorithm. At each time-step, it 
requires the construction of the step operators U{A) 
and V{A) and therefore the diagonalization of h^^ for 
every timc-stcp. As mentioned before, the righthand 
sides of Eqs.fTOl) to ((22l) depend on the Green functions 
at the times T + A which are not known at time T. We 
therefore carry out the time-step T ^ T + A twice. The 
procedure is as follows: 

(1) The collision integrals and h^^ at time T are 
calculated from the Green functions for times t, t' < T. 

(2) A step in the Green function G(T) G{T + A) is 
taken according to Eqs.dMl), dHl), dZZD and 

(3) New collision integrals /f(T + A,t),I^{t,T + 
A), /I {T + A,-iT) a,nd I^{-iT,T+ A) are calculated by 
inserting the new Green functions for times t,t' <T + A 
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into Eqs.dMl) to (HZl). 

(4) The values of the coUision integrals and 
the Hartree-Fock self-energy are approximated 
by / (/(T) + I{T + A))/2 and = 
(E^^(r) + E^^(T + A))/2 where I{T) and I{T + A) 
are the collision terms calculated under points (1) and 
(3). 

(5) The Green function is then propagated from 
G(T) — > G{T + A) using the average values / and 
h^P = h{T + A/2) + in Eqs.dMl), §B, (ED, dZZl) 
and (EH]). 

This concludes the general time-stepping procedure 
for the Green functions. 

We finally consider the calculation of and from 
Eqs. and ([55)1 . As a consequence of the symmetry 
relation ((50l) . we only need to calculate W^{t,t') for 
t < t' . In a time step from T to T -f A we need 
to calculate W<{t,T + A) ior t < T + A from the 
known values of W<{t,T) for t < T. The first term 
on the righthand side of Eg. ((52)) can be calculated 
directly from G<(i,T-h A) and G>(T-h A,t). However, 
the last term X<{t,T + A) of Eq.dSg]) depends on 
the, still undetermined, values W'^{t,T + A). We 
therefore employ an iterative scheme. As a first guess 
for W<{t,T + A) we take W<{t,T + A) = W<{t,T) 
for i < T and W<{T + A,T + A) = W<{T,T). We 
therefore use the values of on the known sides of 
the time square at time T (solid lines in Fig[21) as initial 
guesses for the stepped sides (dotted lines in FigH]) at 
T -|- A. As an initial guess for the value of at the 
new diagonal point (T -I- A, T + A), we take the value at 
the previous diagonal point {T,T). We then calculate 
the quantity X<{t,T + A) for i < T + A and obtain a 
new value for W<{t,T + A) from Eq.^^. This value 



is then inserted back into the righthand side of Eq. (l52)) 
and the process is repeated until convergence is reached. 
Similarly we initiahze W'l (T + A, -ir) = (T, -ir) 
and solve Eq. (|53p in the same manner as for W^. 
This concludes our derivation of the time-stepping al- 
gorithm of the KB equations. The propagation method 
described here has been used in two recent Letters 0, [lH 
where also values for the numerical parameters are 
given. It is clear that the choice of these parameters 
depends strongly on the type of system considered, and 
on the strength of the applied external fields. 



V. SUMMARY AND CONCLUSIONS 

We presented a detailed account of the KB equations 
and discussed in detail their structure, initial conditions 
and symmetries. We developed an algorithm for the 
time-propagation of the KB equations in which the sym- 
metry relations for the Green functions were used to 
reduce the set equations that needed to be solved. In 
two recent Letters ^Z, JLI] we applied the method to the 
case of atoms and molecules in external time-dependent 
fields and to the case of transient transport dynamics of 
double quantum dots. We therefore conclude that time- 
propagation of the KB equations can be used as a practi- 
cal method to calculate the nonequilibrium properties of 
a wide variety of many-body quantum systems, ranging 
from atoms and molecules to quantum dots and quan- 
tum wells. Moreover, the present work can be readily 
extended to other Green function formalisms, such as the 
Nambu formalism [2^, |2^ for superconducting systems. 
Also future extension to bosonic systems is straightfor- 
ward. Work along these lines is in progress. 
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